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I — I 1 Introduction 

< 

A major challenge facing the community working on integral equation based 
solvers for boundary value problems is the construction of efficient discretiza- 
tions of integral operators with singular kernels on curved surfaces. Classic 
approaches such as singularity subtraction, special purpose quadrature, sin- 
gularity cancellation, kernel regularization, and various adaptive strategies 
may work well in many situations but have not yet fully, in three dimensions, 
^ succeeded in unleashing the computational power of integral equation meth- 

yl ods needed for excellence in real- world physics applications [5^ Section 1]. 

Recently, two new promising methods have been launched: the quadrature 
^ by extension (QBX) method which exploits that fields induced by integral 

* operators are often smooth close to the boundaries where their sources are 

^ located [5j and a method relying on a combination of adaptivity, local invert- 

ed ible affine mappings with certain orthogonality properties, and the use of 

precomputed tables of quadrature rules [1] . It seems to be an open question 
what method, or combination of techniques, is best. 
^ This note is about promoting a classic technique for the discretization 

^ of singular integral operators on curved surfaces, namely singularity sub- 

traction. The idea is to use analytical evaluation to a maximum degree and 
split singular (and nearly singular) operators into two parts each - one ill- 
behaved part whose action can be evaluated using high-order analytic prod- 
uct integration, and another more regular part for which purely numerical 
integration is used, compare j2]. Based on this idea we present and imple- 
ment a simple Nystrom scheme for Laplace's equation on tori. Surprisingly 
accurate results are produced. 
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Figure 1: Left: a torus with 6i = 0.5, 62 = 1, pi = 4, p2 = 8 and 32 patches 
Tij. Right: 61 = 0, ^2 = 0.25, pi = 4, p2 = 16 and 64 patches Tij. 

2 Problem formulation and methods 

We consider the interior Dirichlet Laplace problem 

At/(r) = 0, reV, (1) 

U{r)=g{r), rGT, (2) 

where g{r) is a smooth function on the boundary F of a smooth domain V 
in R^. For the solution of (Tj2) we use the double-layer representation 



where is the exterior unit normal of F at position r, da is an element 
of surface area, and is an unknown layer density. An integral equation 
formulation for (T|2) reads 



2.1 Parameterization 

The domain V is taken to be a torus whose surface F is parameterized over 
the square {5 = (51, 52) G M? : — tt < 5i, 52 < vr} as 

r{s) = [g{s) cos(52), g{s) sin(52), S2 sin(5i)] , (5) 

where 

g(s) = 2 + 5i cos(252) + S2 cos(5i) (6) 

and 61 and 62 are shape parameters. The choice 5i = corresponds to the 
standard tori used in ^ Section 3.5]. 

We shall use Nystrom discretization for Q based on composite tensor 
product Gauss-Legendre quadrature. For this, we introduce a sequence of 
mappings pij with z = 1, . . . ,pi and j = 1, . . . ,p2 

Pij{t) = r(^(ti + 2z -pi - l)/pi,7r{t2 + 2j -p2 - l)/p2) • (7) 
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The mapping pij (t) covers a patch Tij of F when mapped from the square 
[t = (ti, ^2) G : -1 < ti,t2 < 1}. The disjoint union of the Tij is T. See 
Figure [1] for two examples. 
Using ([7]) and introducing 

u,j{r,t') = pij(t')-r. (8) 
J«(..0=(^x^)^„.,(r,t'), (9) 

we can rewrite the integral operator in Q as the sum 

PUP2 PUP2 ^ .1 .1 7./^,/X 

E^A,.M^ i: ^/_J_^^.(/)*;d,i. (10) 

2.2 Singularity subtraction 



Nystrom discretization works well for a particular Dij in (10) if r is far 
away from the patch Tij. If r is close to, or on, Tij then the kernel is nearly 
singular, or singular, and something better is needed. Let 



Vi,{r,t') = j2it'k-tk{r))^{t{r)), (11) 
Aij{r,t') = \uijir,t')\^-\vij{r,t')\\ (12) 
where t{r) = p'^^{r). For t' close to t{r) the operator Dij can be expanded 

oo 

Dij = J2Dijk, (13) 

k=0 

See jH Section 1] for a discussion of a similar expansion. 

Our proposed singularity subtraction technique for close to t{r) makes 
use of the split 

Dij^D^ + D°^. (15) 

Here 

K 



with K a small integer and with D^j as in (10). The action of Dj^ is to be 
evaluated using high-order analytic product mtegration and D°j is supposed 
to be sufficiently smooth as to allow for accurate Nystrom discretization. 
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3 Recursive evaluation of integrals 



Computing Dijj^/i{r) in (14) requires the evaluation of an expression of the 
form 



V- r r amntTt2dt[dt', 

(17) 

where a, 6, and c are constants known from (11) and amn are coefficients 



of a polynomial approximating the smooth function Jij{r^t')l\l-{r^t')ji{r'). 
The variable substitution x — at[, xq = ati, y = 6^2, yo = bt2 makes the 



terms in (17) appear as 

1 r «mnX'"i/"dxdi/ 

a^+'b-+^ J -a J-b {{x - xo? + 2c{x - xo){y - yo) + (v - yo?f^^'^ ' 



We now present a scheme for the evaluation of integrals of the form (18). 

Let 

= x2 + 2cxy + y2 (19) 
and define in the Hadamard finite part sense the indefinite integrals 

f f x^y^ dx diy 

Cmnk{x,xo,y,yo,c) ^ / / — VIIZTT^' (20) 

J J \dc[x - xo,y -yo)r+'-/'' 

/x^ dx 
K(x-xo,,-,o)rV2' (21) 

/y^ dy 
u"7 Mfei 1/2 • (22) 
\dc{x -xo,y- yo)r+-^/^ 

Using partial integration, and for given values (x, xq, yo,c)^ one can show 
that when m + n + 1 7^ 2/c holds 

_ mxoC(^_i)^fe + nyoCm(n-i)k + {x - xo)x'^Gnk + {y - yo)y''Fmk 

(23) 

When m + n + 1 — 2k holds 



Cmnk = a:^oC(m-l)nfc + Pki^im - l)C(m_2)„(/c-l) " CnC\m_i)(^n-l)(k-l) 

- x^-'G„(^k-D + c2/"i^(r„-i)(it-i)) , (24) 

Cmnk = yoCm{n-l)k + ((«- " l)C'm(n-2)(fc-l) " CmC(^_i)(„_i)(fc_i) 

- + cx"^G(„_i)(,_i)) , (25) 
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where — ~ c^){2k — 1)). For the evaluation of F^/^ we use 



^00 



mO 



mk 



log {\dc{x -xo,y- yo)|^^^ + (x - Xq) + c{y - yo)^ , 



(x - Xq) + c{y - yo) 
{y-yoy \\dc{x - xo,y -yo)\^-^/'^ 



+ 2(A:-l)Fc 



O(k-l) I 5 



(26) 

> 1, 
(27) 



— (x™ ^|(ic(a; -xo,y - yo)|^^^ - (m - l)(ic(-xo,y - yo)^(m-2)o 
+ (2m-l)(xo-c(y-yo))i^(m-i)o), > 1 , (28) 



1 



(m - 1)F( 



(m-2)(A:-l) 



— X 



m—1 



2k -1 

+ {xo - c{y - yo))F(^rn-i)k , m,k>l. 



\dcix -xo,y- yo 



(29) 



Expression for Gnk obtained by interchanging x ^ y and m n in 
the expressions for F^/c- The recursions for Cmnk^ ^mfc aiid Gnk allow the 
integrals in (17) to be evaluated at a modest computational cost. Note that, 
for this, only Cmnk with m, n > and k > 1 are needed. 



4 Details on the discretization 

We now give precise details on our Nystrom discretization of Q. Aiming at 
10th order convergence we take 10-point composite tensor product Gauss- 
Legendre quadrature (GLIO) as our underlying quadrature scheme. On each 
r^j there will then be a grid of 100 discretization points where the discretized 
density /x is sought. The discretized system Q has 100pip2 unknowns. We 
shall also use a temporary, finer, grid with 256 discretization points on each 
Tij placed according to 16-point composite tensor product Gauss-Legendre 
quadrature (GL16). 



If, for a particular Dij and r in (10), the local parameter t = Pij^ir) is 



such that 3.5 < then the point r is considered far away from Vij and we 
discretize DijiJi{r) using the underlying GLIO scheme. 

If 2 < |t| < 3.5, then r is somewhat close to Vij and we use an extended 
scheme: first /x is interpolated to the finer grid on Vij and then Dijfi(r) 
is discretized using GL16. High-degree polynomial interpolation of smooth 
functions known at Legendre nodes can be very accurate, despite involving 
ill-conditioned Vandermonde systems [3^ Appendix A]. 

If \t\ < 2, then r is close to Tij, the operator Dij is nearly singular or 
weakly singular, and we use the split (15). The discretization is carried out 
on the GL16 grid on F^j, which means that /x has to be interpolated to 
256 points as in the previous paragraph. The operator is discretized 
using GL16. The operator D !^ is discretized using the method of Section [s 
We let m,n = 0, . . . , 15 in |T7| . The 256 coefficients amn are obtainec 
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Convergence with mesh refinement 



5=(0,1) 

5=(0.5,1) 

5=(0,0.25) 
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Number of discretization points 



Figure 2: Relative error in U{r) at the tube center when solving an interior 
Dirichlet Laplace problenn on tori given by ([5|. 



by multiplying the pointwise values of /x at the 256 fine grid points on Tij 
with pointwise values of Jij{r,f)Afj{r,f) and then, in principle, solving a 
Vandermonde system of size 256 x 256. In practice one can obtain the OLrnn 
by solving two 16 x 16 systems with multiple right hand sides. As for the 
optimal number K in ([iS]), it turns out to be related to the polynomial 
degree of the underlying discretization and to the overall mesh refinement 
determined by pi and P2 of ([T]). For small vales pi and p2 and high degree 
quadrature, K should be rather low. We choose K — 1^ that is, we use two 
terms in the sum of (16). 



5 Numerical examples 

Numerical experiments are performed on tori given by ^ using a program 
solely implemented in Matlab and executed on a workstation equipped 
with an IntelXeon E5430 CPU at 2.66 GHz and 32 GB of memory. Three 
different 5 = {81,82) are chosen: 8 = (0,1), 8 = (0.5,1), and 8 = (0,0.25). 
See Figure [1] for illustrations. The boundary condition g{r) in ^ is taken 
as g{r) = l/|r — ri| — l/|r — r2|, with ri = (4, 0, 0) and r2 = (0, 4, 0) for 
8 = (0,1), with n = (4.5,0,0) and = (0,3.5,0) for 8 = (0.5,1), and 
with n = (3.25,0,0) and r2 = (0,3.25,0) for 8 = (0,0.25). The discretized 
system Q is solved iteratively using GMRES. 

Figure [2] shows convergence of U{r), evaluated via a discretization of ([s]), 
at points along the center of the torus tubes. The mesh is refined by in- 
creasing the parameters pi and P2 of Q, keeping p2 — pi for 8 = (0, 1), 
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P2 = 2pi for 6 = (0.5, 1), and p2 = 4pi for 6 — (0, 0.25). A relative residual 
less than 6mach is obtained in between 15 and 20 iterations for reasonably 
resolved systems. The recursion of Section |3] is rather fast. For example, 
with 10,000 discretization points and S = (0, 0.25) only 28 seconds are spent 
doing singularity subtraction. 

One can see in Figure [2] that the initial convergence of U{r) is approx- 
imately 10th order, as expected. As the number of discretization points 
grows, however, the error stemming from the discretization of D°j domi- 
nates and the convergence slows down. Our scheme can, on its own, not 
compete with the mix of techniques presented by Bremer and Gimbutas [Ij. 

6 Conclusion 

This note is about promoting singularity subtraction as a helpful tool in 
the discretization of singular integral operators on curved surfaces. Singular 
and nearly singular kernels are expanded in series whose terms are integrated 
on parametrically rectangular regions using high-order product integration, 
thereby reducing the need for spatial adaptivity and precomputed weights. 
A simple scheme is presented and an application to the interior Dirichlet 
Laplace problem on some tori gives around ten digit accurate results using 
only two expansion terms and a modest programming- and computational 
effort. Further development, including modifications as to allow for para- 
metrically triangular regions, is needed before the technique may find its 
way into competitive solvers. 
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